L3_ Segmentation Methods

Load Data In

load("./data/elk_processed.rda")

Select Individuals (Non-Residential)

library(ggplot2)

ggplot(data=elk_processed, aes(x=lon, y=lat)) +
  geom_path(size=0.5, color="darkgrey") +
  theme_classic() +
  facet_wrap(~id, scale="free", ncol=3)

library(dplyr)

elk_mig <- elk_processed |>
  filter(id %in% c("GP2", "YL25", "YL29", "YL73", "YL77", "YL78"))

ggplot(data=elk_mig, aes(x=datetime, y=lat)) +
  geom_path(size=0.5, color="darkgrey") +
  theme_classic() +
  facet_wrap(~id, scale="free", ncol=2)

# Data Regularization (AniMotum)

install from my R-universe repository

need latest TMB and Matrix package versions also

install.packages("aniMotum", 
                 repos = c("https://cloud.r-project.org",
                           "https://ianjonsen.r-universe.dev"),
                 dependencies = TRUE)


remotes::install_version('TMB', version = '1.9.10')
install.packages("Matrix")

library(aniMotum)

data should be structured with columns: id, date (POSIXct), lc, lon, lat and (optional) additional locational error information smaj, smin, eor (Argos tags)

lc (location class) can include the following values: 3, 2, 1, 0, A, B, Z, G, or GL. The latter two are for GPS locations and ‘Generic Locations’, respectively. Class Z values are assumed to have the same error variances as class B. By default, class G (GPS) locations are assumed to have error variances 10x smaller than Argos class 3 variances, but unlike Argos error variances the GPS variances are the same for longitude and latitude

elk_mig2 <- elk_mig |>
  mutate(lc = "G", date = datetime) |>
  dplyr::select(id, date, lc, lon, lat)

check sampling rate

dtime <- function(t, ...) {difftime(t[-1], t[-length(t)], ...) %>% as.numeric}

sample.rate <- elk_mig2 |> 
  group_by(id) |>
  arrange(date) |>
  mutate(dtime = c(0,round(dtime(date, units ="hours")))) |>
  mutate(mode_dtime = as.numeric(names(sort(table(dtime), decreasing=TRUE))[1])) |>
  data.frame() |>
  ungroup()

barplot(table(sample.rate$dtime))

use sample rate of 2 hours

vmax: m/s

time.step: hrs

random walk: better for large data gaps than crw

elk_ssm_fit <- fit_ssm(elk_mig2,
               vmax = 20,
               model = "rw",
               time.step = 2,
               control = ssm_control(verbose = 0))

summary(elk_ssm_fit)

##  Animal id Model Time n.obs n.filt n.fit n.pred n.rr converged    AICc
##        GP2    rw    2  6082      0  6082   3546    .      TRUE  4542.8
##       YL25    rw    2  6413      1  6412   4664    .      TRUE 14777.5
##       YL29    rw    2  5405      3  5402   3806    .      TRUE 10070.9
##       YL73    rw    2  9614      0  9614   3076    .      TRUE  8472.5
##       YL77    rw    2  2802      0  2802   3230    .      TRUE 13299.3
##       YL78    rw    2  3929      0  3929   2810    .      TRUE  6841.3
## 
## --------------
## GP2 
## --------------
##  Parameter Estimate Std.Err
##      rho_p   0.0622  0.0128
##    sigma_x   0.4546  0.0041
##    sigma_y    0.467  0.0042
##      rho_o  -0.5354  1.7234
##      tau_x   0.0266  0.0176
##      tau_y   0.0288  0.0247
## 
## --------------
## YL25 
## --------------
##  Parameter Estimate Std.Err
##      rho_p    0.219  0.0132
##    sigma_x   0.6061  0.0054
##    sigma_y   0.4665  0.0049
##      rho_o  -0.9189  0.0762
##      tau_x    0.159  0.0317
##      tau_y   0.4218  0.0297
## 
## --------------
## YL29 
## --------------
##  Parameter Estimate Std.Err
##      rho_p   0.1668  0.0133
##    sigma_x   0.5614  0.0054
##    sigma_y   0.5372  0.0052
##      rho_o  -0.8604  0.7478
##      tau_x   0.0531  0.0216
##      tau_y   0.0492  0.0266
## 
## --------------
## YL73 
## --------------
##  Parameter Estimate Std.Err
##      rho_p  -0.0781  0.0101
##    sigma_x   0.6058  0.0044
##    sigma_y   0.5927  0.0043
##          .        .       .
##      tau_x   0.0519  0.0178
##      tau_y   0.0542  0.0226
## 
## --------------
## YL77 
## --------------
##  Parameter Estimate Std.Err
##          .        .       .
##    sigma_x   0.5165  0.0059
##    sigma_y    0.464  0.0051
##          .        .       .
##          .        .       .
##          .        .       .
## 
## --------------
## YL78 
## --------------
##  Parameter Estimate Std.Err
##      rho_p  -0.0461  0.0161
##    sigma_x    0.513  0.0058
##    sigma_y   0.5159  0.0059
##      rho_o  -0.6763  0.7341
##      tau_x   0.0571  0.0243
##      tau_y    0.091  0.0404

extractPredictions <- function(ssm_fit){
  predictions <- data.frame(ssm_fit$predicted) |>
    sf::st_as_sf() |>
    sf::st_transform(4326)
  
  coords <- sf::st_coordinates(predictions)
  
  predictions$lon <- coords[, 1]
  predictions$lat <- coords[, 2]
  
  new_data <- predictions |> 
    data.frame() |>
    dplyr::select(id, date, lon, lat)
  
  return(new_data)
}

elk_ssm_list <- elk_ssm_fit$ssm

elk_regdata_list <- lapply(elk_ssm_list,
                      extractPredictions)

library(ggplot2)

elk_id <- split(elk_mig2, elk_mig2$id)

track_plots <- c()
for(i in c(1:length(elk_id))){
  track_plots[[i]] <- ggplot()+
    geom_path(data = elk_regdata_list[[i]], aes(x = lon, y = lat), size = 0.7, color = "red")+
    geom_point(data = elk_id[[i]], aes(x = lon, y = lat), size=1, color="black")+
    theme_classic()+
    xlab("Longitude")+ylab("Latitude")+
    facet_wrap(~id,scales="free")
}

track_plots

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

elk_regdata <- do.call("rbind", elk_regdata_list)

Calculate Movement Metrics

library(sf)

elk_reg_utm <- elk_regdata |>
  st_as_sf(coords = c("lon","lat"), crs = 4326) |>
  st_transform(32611)

coords <- st_coordinates(elk_reg_utm)

elk_regdata$X <- coords[,1]
elk_regdata$Y <- coords[,2]
elk_regdata$Z <- elk_regdata$X + 1i*elk_regdata$Y

getMovementMetrics <- function(dataframe){
  
  move_step <- diff(dataframe$Z)  

  time_step <- as.numeric(difftime(dataframe$date[-1], dataframe$date[-length(dataframe$date)], "days"))

  absolute_turnangle <- Arg(move_step)

  relative_turnangle <- diff(absolute_turnangle)

  step_length_km <- Mod(move_step)/1000
  
  dataframe$time_step_days <- c(NA, time_step)

  dataframe$step_length_km <- c(NA, step_length_km)

  dataframe$relative_turnangle <- c(NA, NA, relative_turnangle)
  
  return(dataframe)
}

elk_reg_id <- split(elk_regdata, elk_regdata$id)

elk_reg_id2 <- lapply(elk_reg_id,
                            getMovementMetrics)

Lavielle

The Lavielle’s method identifies breaking points in a timeseries by using a moving window through the timeseries and a penalized contrast function which determines the optimal number of segments in the timeseries by minimizing the penalized contrast function (Lavielle 2005)

The method examine the variation in the mean and/or variance of a timeseries and identify the optimal number of break points, by fitting the mean, variance or both to the timeseries depending on the number of segments (K). As K increases, the fit increases. However, there should be a number of segments K after which adding more segments does not contributes significantly to improving the fit. The penalized contrast function adjust the trade-off between fit and K

library(adehabitatLT)

str(elk_reg_id2[["YL73"]])

## 'data.frame':    3076 obs. of  10 variables:
##  $ id                : chr  "YL73" "YL73" "YL73" "YL73" ...
##  $ date              : POSIXct, format: "2004-02-20 00:00:00" "2004-02-20 02:00:00" ...
##  $ lon               : num  -116 -116 -116 -116 -116 ...
##  $ lat               : num  51.8 51.7 51.7 51.7 51.7 ...
##  $ X                 : num  602201 599949 599230 599220 599239 ...
##  $ Y                 : num  5734480 5733904 5733606 5733601 5733364 ...
##  $ Z                 : cplx  6e+05+5.73e+06i 6e+05+5.73e+06i 6e+05+5.73e+06i ...
##  $ time_step_days    : num  NA 2 2 2 2 2 2 2 2 2 ...
##  $ step_length_km    : num  NA 2.3245 0.778 0.0113 0.2374 ...
##  $ relative_turnangle: num  NA NA 0.1424 0.0303 1.2283 ...

ggplot(data = elk_reg_id2[["YL73"]], aes(x = date, y = lat))+
  geom_path(color="darkgrey") +
  geom_point(size=0.5) +
  theme_classic() +
  xlab("Date") + ylab("Latitude")

decide a minimum number of relocations within a segment (Lmin), the maximum number of segments we want the function to try (Kmax) and if it should look at the difference in mean, variance, or both (mean, var, meanvar)

elk_example_data <- elk_reg_id2[["YL73"]]

elk_lavielle_example <- lavielle(elk_example_data$lat, Lmin = 50, Kmax = 10, type = "mean")

chooseseg(elk_lavielle_example)

##   K            D         Jk
## 1 1          Inf 3075.00000
## 2 2  4.914028021  860.63701
## 3 3  1.313698068  282.70349
## 4 4  0.379474355  142.24697
## 5 5  0.000609291  128.15990
## 6 6  0.001766675  114.27573
## 7 7  0.021777072  100.97988
## 8 8 -0.014404272   94.93606
## 9 9  0.013929530   84.09544

elk_lavielle_example_path <- findpath(elk_lavielle_example, 3)

breakpoints <- unlist(elk_lavielle_example_path)

breakpoints

## [1]    1  892  893 2924 2925 3076

breakpoints_df <- data.frame(Row_ID = breakpoints, segment = c(1,1,2,2,3,3))

elk_example_data$Row_ID <- cumsum(rep(1, nrow(elk_example_data)))

elk_example_data_lav <- merge(elk_example_data, breakpoints_df, by="Row_ID", all.x=TRUE)

elk_example_data_lav <- tidyr::fill(elk_example_data_lav, segment, .direction = 'down')

head(elk_example_data_lav)

##   Row_ID   id                date       lon      lat        X       Y
## 1      1 YL73 2004-02-20 00:00:00 -115.5194 51.75189 602200.8 5734480
## 2      2 YL73 2004-02-20 02:00:00 -115.5522 51.74711 599949.0 5733904
## 3      3 YL73 2004-02-20 04:00:00 -115.5627 51.74456 599230.3 5733606
## 4      4 YL73 2004-02-20 06:00:00 -115.5628 51.74452 599220.1 5733601
## 5      5 YL73 2004-02-20 08:00:00 -115.5626 51.74239 599239.2 5733364
## 6      6 YL73 2004-02-20 10:00:00 -115.5575 51.73933 599602.2 5733031
##                 Z time_step_days step_length_km relative_turnangle segment
## 1 602201+5734480i             NA             NA                 NA       1
## 2 599949+5733904i              2     2.32454229                 NA       1
## 3 599230+5733606i              2     0.77796897         0.14237354       1
## 4 599220+5733601i              2     0.01127217         0.03027534       1
## 5 599239+5733364i              2     0.23737279         1.22829486       1
## 6 599602+5733031i              2     0.49274756         0.74708438       1

First Passage Time

The first passage time (FPT) is a parameter often used to describe the scale at which patterns occur in a trajectory. For a given scale r, it is defined as the time required by the animals to pass through a circle of radius r.

Fauchald & Tveraa (2003) proposed another use of the FPT. Instead of computing the mean of FPT, they propose the use of the variance of the log(FPT). This variance should be high for scales at which patterns occur in the trajectory (e.g. area restricted search). This method is often used to determine the scale at which an animal seaches for food

library(mapview)

elk_example_track <- elk_example_data |>
  st_as_sf(coords=c("lon","lat"), crs= 4326) |>
  st_union() |>
  st_cast("LINESTRING") 

mapview(elk_example_track)

transform the data into ltraj object

elk_example_sp <- elk_example_data |>
  st_as_sf(coords=c("lon","lat"), crs= 4326) |> 
  st_transform(32611) |>
  st_cast("POINT") %>% as_Spatial()

elk_example_traj <- as.ltraj(coordinates(elk_example_sp), date=elk_example_data$date, id = elk_example_data$id)

summary(elk_example_traj)

##     id burst nb.reloc NAs date.begin            date.end
## 1 YL73  YL73     3076   0 2004-02-20 2004-11-02 06:00:00

head(elk_example_traj[[1]])

##                  x       y                date          dx          dy
## YL73.X    602200.8 5734480 2004-02-20 00:00:00 -2251.83611 -576.828406
## YL73.X.3  599949.0 5733904 2004-02-20 02:00:00  -718.61809 -298.033141
## YL73.X.6  599230.3 5733606 2004-02-20 04:00:00   -10.27673   -4.631476
## YL73.X.9  599220.1 5733601 2004-02-20 06:00:00    19.18594 -236.596153
## YL73.X.12 599239.2 5733364 2004-02-20 08:00:00   362.94768 -333.270359
## YL73.X.15 599602.2 5733031 2004-02-20 10:00:00   116.54171    8.977010
##                 dist   dt      R2n  abs.angle  rel.angle
## YL73.X    2324.54229 7200        0 -2.8908256         NA
## YL73.X.3   777.96897 7200  5403497 -2.7484521 0.14237354
## YL73.X.6    11.27217 7200  9588981 -2.7181767 0.03027534
## YL73.X.9   237.37279 7200  9658265 -1.4898819 1.22829486
## YL73.X.12  492.74756 7200 10016404 -0.7427975 0.74708438
## YL73.X.15  116.88694 7200  8853351  0.0768765 0.81967397

plot(elk_example_traj)

radii - numeric vector for the radii of the circles

units - time units of the results

elk_example_fpt <- fpt(elk_example_traj, radii = seq(50, 10000), units="hours")

peak in variance at spatial scale of interest

varlogfpt(elk_example_fpt, graph=TRUE)

elk_example_fpt <- fpt(elk_example_traj, radii = 10000, units="hours")

plot(elk_example_fpt, scale = 9000, warn = FALSE)

elk_example_data$FPT <- elk_example_fpt[[1]]$r1

ggplot() +
  geom_path(data = elk_example_data, aes(x = date, y = scale(FPT)), color="darkgrey", size = 1) +
  geom_path(data = elk_example_data, aes(x = date, y = scale(lat)), color ="orange", alpha = 0.6, size = 1) +
  geom_point(data = elk_example_data, aes(x = date, y = scale(FPT))) +
  theme_classic() +
  ylab("Scaled FPT (grey) vs Latitude (orange)")

Behavioral Change Point Analysis

The “smoove” package allows for multiple continuous-time movement models (correlated velocity models or CVMs) to be fit to irregular tracking data to identify “modes” of behavior. Each model comes with it’s own set of parameters and best models are determined using AIC/BIC and maximum likelihood. CVM model options include: Uncorrelated CVM (UCVM), Advective CVM (AVCM), Rotational CVM (RCVM), and Rotational-Advective CVM (RACVM). A single-change point can also be identified assuming just a UCVM process and sudden changes in the UCVM parameters, time scale, and root mean square speeds.

library(devtools)

install_github("EliGurarie/smoove")

library(smoove)

elk_utm <- elk_regdata |>
  st_as_sf(coords = c("lon","lat"), crs = 4326) |>
  st_transform(32611)

coords <- st_coordinates(elk_utm)

elk_regdata$X <- coords[,1]

elk_regdata$Y <- coords[,2]

elk_regdata$Z <- elk_regdata$X + 1i*elk_regdata$Y

elk_example <- elk_regdata |>
  dplyr::select(id, date, X, Y, Z) |>
  dplyr::filter(id=="YL73")

with(elk_example, scan_track(x = X, y = Y, main = "YL73"))

specify window size (how long does behavior last?)

try: 200 days

takes awhile to run

elk_sweep <- with(elk_example, sweepRACVM(Z = Z, T = date, windowsize = 200, time.unit = "days", windowstep = 30, model = "RACVM", progress = TRUE))

parallel example

library(doParallel)

cl <- parallel::makeCluster(parallel::detectCores())

doParallel::registerDoParallel(cl)

elk_sweep <- with(elk_example, sweepRACVM(Z = Z, T = date, windowsize = 200, time.unit = "days", windowstep = 30, model = "RACVM", progress = TRUE, .parallel = TRUE))

colors represent the relative log-likelihood profile for a single window

looking for distinct peaks (significant change points)

plotWindowSweep(elk_sweep, main = "YL73")

Can specify the clusterwidth here- I used widths ranging from 1-6, based on warnings given with fitting the function for each id

I used all models as candidates (UCVM, ACVM, RCVM, RACVM)

I additionally used AIC to determine models, as BIC can be overly conservative

elk_cp <- findCandidateChangePoints(windowsweep = elk_sweep, clusterwidth = 6) |>
  getCPtable(modelset = "all",criterion="AIC")

estimate phases

elk_phases <- elk_cp |>
  estimatePhases()

##   phase               start                 end model       eta tau       rms
## 1     I 2004-02-20 00:00:00 2004-05-05 10:00:00  UCVM 10747.467   0 10747.467
## 2    II 2004-05-05 10:00:00 2004-11-02 06:00:00  UCVM  8627.125   0  8627.125

elk_phases_summ <- summarizePhases(elk_phases)

str(elk_phases_summ)

## 'data.frame':    2 obs. of  7 variables:
##  $ phase: chr  "I" "II"
##  $ start: POSIXct, format: "2004-02-20 00:00:00" "2004-05-05 10:00:00"
##  $ end  : POSIXct, format: "2004-05-05 10:00:00" "2004-11-02 06:00:00"
##  $ model: chr  "UCVM" "UCVM"
##  $ eta  : num  10747 8627
##  $ tau  : num  0 0
##  $ rms  : num  10747 8627

elk_phases_data <- merge(elk_example, elk_phases_summ, all = TRUE, by.x = c("date"), by.y = c("end")) |> arrange(date) |>
  tidyr::fill(phase, start, model, eta, tau, rms, .direction = "updown")

str(elk_phases_data)

## 'data.frame':    3077 obs. of  11 variables:
##  $ date : POSIXct, format: "2004-02-20 00:00:00" "2004-02-20 02:00:00" ...
##  $ id   : chr  "YL73" "YL73" "YL73" "YL73" ...
##  $ X    : num  602201 599949 599230 599220 599239 ...
##  $ Y    : num  5734480 5733904 5733606 5733601 5733364 ...
##  $ Z    : cplx  6e+05+5.73e+06i 6e+05+5.73e+06i 6e+05+5.73e+06i ...
##  $ phase: chr  "I" "I" "I" "I" ...
##  $ start: POSIXct, format: "2004-02-20 00:00:00" "2004-02-20 00:00:00" ...
##  $ model: chr  "UCVM" "UCVM" "UCVM" "UCVM" ...
##  $ eta  : num  10747 10747 10747 10747 10747 ...
##  $ tau  : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ rms  : num  10747 10747 10747 10747 10747 ...

ggplot(data = elk_phases_data, aes(x = date, y = Y, color = model))+
    geom_path(size=1.5)+
    xlab("Datetime")+ylab("Latitude")+
    theme_classic()

ggplot(data = elk_phases_data, aes(x = date, y = Y, color = phase))+
    geom_path(size=1.5)+
    xlab("Datetime")+ylab("Latitude")+
    theme_classic()

Hidden Markov Model

library(momentuHMM)

elk_hmm_prep <- prepData(elk_example, type = "UTM",
                                    coordNames = c("X","Y")) 

head(elk_hmm_prep)

##        ID       step      angle   id                date               Z
## 1 Animal1 2324.54229         NA YL73 2004-02-20 00:00:00 602201+5734480i
## 2 Animal1  777.96897 0.14237354 YL73 2004-02-20 02:00:00 599949+5733904i
## 3 Animal1   11.27217 0.03027534 YL73 2004-02-20 04:00:00 599230+5733606i
## 4 Animal1  237.37279 1.22829486 YL73 2004-02-20 06:00:00 599220+5733601i
## 5 Animal1  492.74756 0.74708438 YL73 2004-02-20 08:00:00 599239+5733364i
## 6 Animal1  116.88694 0.81967397 YL73 2004-02-20 10:00:00 599602+5733031i
##          x       y
## 1 602200.8 5734480
## 2 599949.0 5733904
## 3 599230.3 5733606
## 4 599220.1 5733601
## 5 599239.2 5733364
## 6 599602.2 5733031

set any step lengths with a value of 0 to 1 (can’t have 0’s)

elk_hmm_prep$step[elk_hmm_prep$step == 0] <- 1

define priors and prior distributions

1 - state 1

2 - state 2

stepMean0 <- c(m1 = 100, m2 = 4000)
stepSD0 <- c(sd1 = 50, sd2 = 1000)
angleCon0 <- c(rho1  = 0.1, rho2 = 0.8)

stateNames <- c("resident","transit")

set distributions for variables (step lengths and turning angles)

set priors (mean and sd for steps, concentration for angle)

see ?dgamma and ?dwrpcauchy

dist <- list(step = "gamma", angle = "wrpcauchy")
Par0 <- list(step=c(stepMean0, stepSD0), angle = c(angleCon0))

nbStates = number of states to identify

elk_hmm_fit <- fitHMM(data = elk_hmm_prep, nbStates = 2, dist = dist, 
                      Par0 = Par0, stateNames = stateNames)

print(elk_hmm_fit)

## Value of the maximum log-likelihood: -27262.61 
## 
## 
## step parameters:
## ----------------
##      resident  transit
## mean 215.5848 902.3420
## sd   219.7929 767.8621
## 
## angle parameters:
## -----------------
##                    resident   transit
## mean           0.000000e+00 0.0000000
## concentration 1.458213e-172 0.4219822
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##                1 -> 2     2 -> 1
## (Intercept) -1.338506 -0.5050826
## 
## Transition probability matrix:
## ------------------------------
##           resident   transit
## resident 0.7922442 0.2077558
## transit  0.3763470 0.6236530
## 
## Initial distribution:
## ---------------------
##     resident      transit 
## 6.537404e-06 9.999935e-01

viterbi formula to decode the transition probabilities for most likely states at each time point/location

hmm_states <- viterbi(elk_hmm_fit)

str(hmm_states)

##  int [1:3076] 2 2 1 1 1 1 1 2 2 1 ...

what does this code do?

M <- elk_hmm_fit$mle$gamma

step <- elk_hmm_fit$mle$step

angle <- elk_hmm_fit$mle$angle[2,] 

what does this code do? proportion of locs in each state?

probs <- c(s2 = M[1,2]/(M[1,2] + M[2,1]), 
              s1 = M[2,1]/(M[1,2] + M[2,1]))

probs

##        s2        s1 
## 0.3556837 0.6443163

elk_example$state <- hmm_states

layout(cbind(c(1,1),2:3))
par(bty = "l", mar = c(2,2,2,2))

with(elk_example, {
  plot(X, Y, asp =1, col = c("orange","blue")[state], pch = 19, cex = 0.7)
  segments(X[-length(X)], Y[-length(Y)], 
           X[-1], Y[-1], col = c("orange","blue")[state[-length(state)]])
  plot(date, X, col = c("orange","blue")[state], pch = 19, cex = 0.7)
  segments(date[-length(X)], Y[-length(Y)], 
           date[-1], Y[-1], col = c("orange","blue")[state[-length(state)]])
  plot(date, Y, col = c("orange","blue")[state], pch = 19, cex = 0.7)
  segments(date[-length(X)], Y[-length(Y)], 
           date[-1], Y[-1], col = c("orange","blue")[state[-length(state)]])
})

elk_example_sf <- elk_example |>
  st_as_sf(coords=c("X","Y"), crs= 32611) |>
  st_transform(4326) |>
  mutate(state = as.character(state))

elk_example_track <- elk_example_sf |>
  summarize(do_union=FALSE) |> 
  st_cast("LINESTRING")

mapview(elk_example_track, color="darkgrey") +
  mapview(elk_example_sf, zcol="state", col.regions=c("orange","blue"))

refit and add covariate (latitude)

formula <- ~y

elk_hmm_fit2 <- fitHMM(data = elk_hmm_prep, nbStates = 2, dist = dist, 
                      Par0 = Par0, stateNames = stateNames, formula = formula)

print(elk_hmm_fit2)

## Value of the maximum log-likelihood: -38803.7 
## 
## 
## step parameters:
## ----------------
##      resident transit
## mean      100    4000
## sd         50    1000
## 
## angle parameters:
## -----------------
##               resident transit
## mean               0.0     0.0
## concentration      0.1     0.8
## 
## Regression coeffs for the transition probabilities:
## ---------------------------------------------------
##             1 -> 2 2 -> 1
## (Intercept)   -1.5   -1.5
## y              0.0    0.0
## 
## Transition probability matrix (based on mean covariate values):
## ---------------------------------------------------------------
##           resident   transit
## resident 0.8175745 0.1824255
## transit  0.1824255 0.8175745
## 
## Initial distribution:
## ---------------------
## resident  transit 
##      0.5      0.5

hmm_states <- viterbi(elk_hmm_fit2)

elk_example$state <- hmm_states

layout(cbind(c(1,1),2:3))
par(bty = "l", mar = c(2,2,2,2))

with(elk_example, {
  plot(X, Y, asp =1, col = c("orange","blue")[state], pch = 19, cex = 0.7)
  segments(X[-length(X)], Y[-length(Y)], 
           X[-1], Y[-1], col = c("orange","blue")[state[-length(state)]])
  plot(date, X, col = c("orange","blue")[state], pch = 19, cex = 0.7)
  segments(date[-length(X)], Y[-length(Y)], 
           date[-1], Y[-1], col = c("orange","blue")[state[-length(state)]])
  plot(date, Y, col = c("orange","blue")[state], pch = 19, cex = 0.7)
  segments(date[-length(X)], Y[-length(Y)], 
           date[-1], Y[-1], col = c("orange","blue")[state[-length(state)]])
})

add a third state

stepMean0 <- c(m1 = 100, m2 = 4000, m3 = 1000)
stepSD0 <- c(sd1 = 50, sd2 = 1000, sd3 = 500)
angleCon0 <- c(rho1  = 0.1, rho2 = 0.8, rho3 = 0.5)

stateNames <- c("resident","transit", "other")

dist <- list(step = "gamma", angle = "wrpcauchy")
Par0 <- list(step=c(stepMean0, stepSD0), angle = c(angleCon0))

formula <- ~y

elk_hmm_fit3 <- fitHMM(data = elk_hmm_prep, nbStates = 3, dist = dist, 
                      Par0 = Par0, stateNames = stateNames, formula = formula)

hmm_states <- viterbi(elk_hmm_fit3)

elk_example$state <- hmm_states

layout(cbind(c(1,1),2:3))
par(bty = "l", mar = c(2,2,2,2))

with(elk_example, {
  plot(X, Y, asp =1, col = c("orange","blue","green")[state], pch = 19, cex = 0.7)
  segments(X[-length(X)], Y[-length(Y)], 
           X[-1], Y[-1], col = c("orange","blue", "green")[state[-length(state)]])
  plot(date, X, col = c("orange","blue", "green")[state], pch = 19, cex = 0.7)
  segments(date[-length(X)], Y[-length(Y)], 
           date[-1], Y[-1], col = c("orange","blue", "green")[state[-length(state)]])
  plot(date, Y, col = c("orange","blue", "green")[state], pch = 19, cex = 0.7)
  segments(date[-length(X)], Y[-length(Y)], 
           date[-1], Y[-1], col = c("orange","blue", "green")[state[-length(state)]])
})

elk_example_sf <- elk_example |>
  st_as_sf(coords=c("X","Y"), crs= 32611) |>
  st_transform(4326) |>
  mutate(state = as.character(state))

elk_example_track <- elk_example_sf |>
  summarize(do_union=FALSE) |> 
  st_cast("LINESTRING")

mapview(elk_example_track, color="darkgrey") +
  mapview(elk_example_sf, zcol="state", col.regions=c("orange","blue", "green"))